Numerical method for binary black hole/neutron star initial data: Code test 
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A new numerical method to construct binary black hole/neutron star initial data is presented. The 
method uses three spherical coordinate patches; Two of these are centered at the binary compact 
objects and cover a neighborhood of each object; the third patch extends to the asymptotic region. 
As in the Komatsu-Eriguchi-Hachisu method, nonlinear elliptic field equations are decomposed into 
a flat space Laplacian and a remaining nonlinear expression that serves in each iteration as an 
effective source. The equations are solved iteratively, integrating a Green's function against the 
effective source at each iteration. Detailed convergence tests for the essential part of the code 
are performed for a few types of selected Green's functions to treat different boundary conditions. 
Numerical computation of the gravitational potential of a fluid source, and a toy model for a 
binary black hole field are carefully calibrated with the analytic solutions to examine accuracy and 
convergence of the new code. As an example of the application of the code, an initial data set for 
binary black holes in the fsenberg-Wilson-Mathews formulation is presented, in which the apparent 
horizons are located using a method described in Appendix A. 



I. INTRODUCTION 

Inspiral and merger simulations to produce accurate 
gravitational waveforms are essential for constructing 
waveform templates for analysis of data from laser in- 
terferometric detectors. The ground based interferome- 
ters, such as Advanced LIGO or LCGT may detect grav- 
itational waves from the inspiral of M ~ 10A/q binary 
black holes within z ~ 4, while space based interferomet- 
ric detectors such as LISA or DECIGO may detect the 
inspiral of lO 6 Af supermassive binary black holes, and 
may discover intermediate mass binary black holes with 
M ~ 1O 3 M . 

Initial data set for binary black holes with a variety 
of black hole parameters, such as binary mass ratio, and 
black hole spins, or the binary black hole-neutron star 
data will become more useful considering the remarkable 
progress made recently for the inspiraling binary black 
holes simulations up to a few orbits near the innermost 
stable circular orbits 0, 0, S HIE IE 0, H] • Several groups 
have been achieved to construct binary black hole initial 
data successfully d EpJLl El 0, 03, El El El , 
(for earlier works, see |19j'). 

In this paper, we introduce a new numerical method 
suitable for computing accurate initial data sets for bi- 
nary black holes, black hole-neutron star binaries, and 
binary neutron star systems. Initial data sets of these 
kinds are calculated from the Einstein equation written 
in the form of nonlinear elliptic equations for metric com- 
ponents. Each equation can be written as a Poisson equa- 
tion with a nonlinear source. Our new Poisson solver 
is patterned after the Komatsu-Eriguchi-Hachisu (KEH) 
method [201 ] . widely used to compute rotating neutron 
stars and more recently to compute, binary neutron stars 
[2l| . The KEH method uses Green's formula to write the 
field equations in equivalent integral forms and iteratively 
solve them using spherical coordinates and angular har- 
monics. The set of equations is discretized by a standard 



finite difference scheme. 

To extend the method to handle our wider class of 
binary configurations, we introduce three spherical co- 
ordinate patches. Two are centered at each hole and 
extend outward to a finite radius, larger than the gravi- 
tational radius but small enough that the two coordinate 
patches do not overlap. A third coordinate patch, cov- 
ering the rest of a spacelike hypersurface, extends to the 
the asymptotic region and overlaps each of the other two 
patches. There are two important features of our new 
code: (1) The number of multipolcs in the coordinate 
patch centered at the orbital center can be reduced to 
I < 10 since the sizes of patches for compact objects is 
extended to about a half of the binary separation. (2) 
The data between those patches are communicated only 
at the boundary of those patches to minimize the amount 
of data to interpolate from one to the other. These novel 
features result in an efficient code that retains high res- 
olution even near the compact objects where the field is 
strong. The method has two additional significant ad- 
vantages: Coding is relatively simple, and the iteration 
converges robustly. 

This paper is organized as follows. In section [Til after 
briefly reviewing the initial value formulation, we intro- 
duce our choice of coordinate systems and the formu- 
lation of our Poisson solver. Formulas for multipolc ex- 
pansion of Green's function used in the Poisson solver are 
described in Appendix [Bl In section HTT1 we describe our 
numerical methods, including finite differencing and our 
iteration procedure. In the relating Appendices [D] and 
lEl the convergence of iteration is discussed further. The 
results of detailed convergence tests are presented in IIV1 
and an example of binary black hole initial data is dis- 
played in [V] The concrete form of the nonlinear elliptic 
equations for the initial value problem is given in Ap- 
pendix [C] and a method to locate the apparent horizons 
in the initial data is described in Appendix [Al 
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II. METHOD FOR BINARY BLACK 
HOLE/NEUTRON STAR INITIAL DATA 

Our new numerical method is applicable for various 
formulations including spatially confomal flat initial data 
22 1, the Isenberg- Wilson-Mathews (IWM) formulation 
23 , [H, [2|| , and waveless approximation [26] . We intro- 
duce the IWM formulation used for a test calculation of 
our new code. In this formulation, four constraints and 
the spatial trace of the Einstein equation are solved, af- 
ter choosing the trace of the extrinsic curvature and the 
conformal three metric. We then explain the choice of 
coordinates and the form of Green's functions used in 
our version of the KEH method. 



A. Formulation 

We consider a globally hyperbolic spacctime A4, fo- 
liated by a family of spacelike hypersurfaces £4. The 
binary black hole initial data is constructed on a slice 
E = So. The unit future-pointing normal to £4 will be 
denoted by n a — —aV a t and the metric is written in 
3+1 form, 

ds 2 — g^dx^dx" 

= -a 2 dt 2 + 7lJ (dx i + /3 i dt){dx j + /3 j dt), (1) 

in a chart {t,x 1 }, where 7ab(i) is the 3-metric on £ t , a 
the lapse and (3 a the shift. The 3-metric 7 a t, is induced 
by the projection tensor to the hypersurfaces £ t 

lap = 9ap + n a np. (2) 

The extrinsic curvature of the foliations is defined by 



K al 3 = --£ n 7a0, 



(3) 



where K a p satisfies K a pn a = 0. With the spatial indices, 
the spatial tensor K a p is written 



Kab = -— {dtlab - £/37q&) • 
la 



(4) 



Denoting the tracefree part of K ab by A ab and its trace 
by K := K a a , we have 



and 



A a b = Kab - -^JabK, 



A a b = — ^ f £n7afc — g labl^ £njcd 



(5) 



(6) 



£ n operating to the spatial metric is understood as £„ = 
— (dt — £/3), where the £p is the Lie derivative defined on 
£. 



Projecting one index of the Einstein equation normal 
to the hypersurface £, G a pn a — 0, yields the Hamilto- 
nian and momentum constraint equations 

2(G Q(3 - ^T a0 )n a nP 

= R- K ab K ab + K 2 - 16^ph = 0, (7) 



(G al3 - 8nT al3 )^ a n a 

= -D b (K ab - j ab K ) + 8iTj a = 



(8) 



where D a is the covariant derivative on E associated with 
the three metric i a b- To satisfy the Hamiltonian con- 
straint on a slice E, we introduce a conformal decom- 
position of the spatial metric, 7 a fc = ip 7 a i>, and solve 
the Hamiltonian constraint for the conformal factor ?/>. 
This prescription leaves the conformal three metric j a b 
unspecified. 

Separating out the trace K, and substituting Eqs. ([5]) 
and ([(j]) in the momentum constraint results in an elliptic 
equation for the shift j3 a . The trace, K, remains unspec- 
ified. The spatial trace of the Einstein equation, 

(G Q/3 - 87^)7^ 

= 2£ n K -\{R + K 2 + 3K ab K ab ) + -D a D a a - 8tt5 
2 a 



0, 



(9) 



can be written as an elliptic equation for the lapse a, 
once one restricts dtK (e.g., by setting to zero dtK or 
the derivative of K along a helical Killing vector) . 

In Appendix[C] we show the explicit form of the elliptic 
equations with nonlinear source for the constraints and 
the spatial trace of the Einstein equation. 



B. Inversion of the Laplacian: Poisson solver 

In each component of the field equations a second-order 
elliptic operator acts on one metric potential. By sepa- 
rating out a flat Laplacian, we write each field-equation 
component in the form 



V 2 $ = 5, 



(10) 



where $ represents a metric potential on a slice E. The 
effective source S involves second derivatives of the met- 
ric, but a convergent iteration is possible because they 
occur in expressions that are o(r~ 3 ) near spatial infinity. 
The flat Laplacian V 2 is separated in spherical coordi- 
nates and inverted by a Poisson solver, and the nonlinear 
equation (fTOf is solved iteratively. Our choice of the Pois- 
son solver is to use the Green's formula, an integral form 
of Eq. (fTU|) . Using the Green's function of the Laplacian 



V 2 G(x,x r ) = ~4tt5(x - x'), 



(11) 



where x and x' are positions, i,i'£FCS, the Green's 
formula is written 



$(a:) = ( G(x,x r )S(x')d 3 x' 

47T Jy 
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+_L J [ G (x,x')V a <S>(x') 
47r JdV 

-<S>{x')V' a G{x,x')]dS' a . (12) 

This formula is valid for any connected space V as long as 
each term is integrable. The function G(x, x') is a sum of 
a Green's function without boundary and a homogeneous 
solution F(x, x') to the Laplace equation, 

G(x,x') = r ^— + F(x,x'), (13) 
| X — x I 

which satisfy 

V 2 i — - — tt = -AttS(x-x'), (14) 
| a; — x'\ 

V 2 F(x,x') = 0. (15) 

Eq. (fT2|) is a formal solution to the Poisson equation fTO]) 
for any G(x, x') that satisfies Eq. (fTT|) . even if the source 
S depends on the field $ nonlincarly. Thereby, the ellip- 
tic equation with a nonlinear source can be solved iter- 
atively using Eq. (fT2|) . We call this iteration the KEH 
iteration hereafter. 

With the surface term included, Eq. (Tj"2")) is an iden- 
tity, valid for any choice of Green's function. Requir- 
ing convergence of the KEH iteration, however, im- 
poses the following key restrictions on that choice: For 
solving a Dirichlet problem, no multipole component of 
X7' a G(x 7 x') can vanish on the entire boundary; similarly, 
for solving a Neumann problem, no multipole compo- 
nent of G(x,x') can vanish on the entire boundary. As 
long as the Green's function satisfies this restriction, it is 
not necessary to construct F(x, x 1 ) appropriate for each 
boundary condition; for example, the Green's function 
without boundary term 1/ \x — x'\ can be used for the 
Neumann problem. 

A reason to use the Poisson solver Eq. (|T2"]) is tied to 
the facts that the spherical coordinates (r, 0, (f>) are suit- 
able for constructing numerical domains for binary black 
holes and neutron stars, and that one can use a multipole 
expansion of the Green's function in these coordinates. 
The Poisson solver turns out to be simple for coding, 
CPU inexpensive and accurate as shown in Sec lIVI In 
the rest of this section, we introduce numerical domains 
for the binary black holes, and discuss the choice of the 
Green's function. 



C. Construction of the computational domain 

To describe the binary black hole/neutron star data, 
we introduce three spherical domains. Two small do- 
mains are centered at the two compact objects surround- 
ing them, while the third domain partly covers the first 
two and extends to the asymptotic region. The large do- 
main is not necessarily positioned at the center of mass 
of the two compact objects. These domains are shown 
schematically in Fig. [TJ In this section, we consider the 




FIG. 1: The computational domain. One central grid and 
two black hole grids with excised regions. 



binary black hole case as an example, and refer to the 
domains around the compact objects as the black hole 
coordinate system (BHCS) and the third one as the cen- 
tral coordinate system (CCS). 

CCS extends to the asymptotic region; practically the 
radius r of the sphere S is set large enough that the 
multipoles of order r~ 2 and higher are negligible. It ex- 
cludes the interiors of the two spheres I\ and I2, which 
are centered at each black hole, and whose radii are taken 
larger than the gravitational radius of each hole but not 
as large as to intersect each other. Therefore, in the do- 
main of CCS, Eq. (|12[) involves (1) the surface integrals 
over a large sphere S where the asymptotic condition of 
each field variable is imposed, (2) the interfaces 1\ and I2 , 
and (3) the volume integral of the source term S between 
these three spheres. 

The first black hole computational domain BHCS-1 ex- 
tends from a sphere to S 0l , both centered at the black 
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hole, and the second one BHCS-2 from Si 2 to S 02 . The 
region inside and Si 2 is excised from the computa- 
tional domain. The code allows the option of dispensing 
with the inner boundary, for a black hole; and for a neu- 
tron star the inner boundary is never used (in this case 
in the BHCS the minimum value of r is zero). Note that 
the spheres I\, Si t , and S 0l are concentric, and the same 
for I 2 , S i2 , and S D2 . 

The boundaries I\, 1%, S 0l , and 5 Q2 are introduced to 
reduce the number of terms in the Legendre expansion of 
the Green's function (|B1[) in CCS. Taking the radii of I\ 
and I2 large enough, the contribution of higher multipoles 
in CCS is included in the surface integrals over I\ and I2 ■ 
Because of this, a small number of multipoles, typically 
£ < 10, is enough to resolve the volume integral of CCS. 
Thus by increasing the radial resolution in BHCS where 
the metric potentials may vary rapidly, we can compute 
an accurate solution without having a high resolution in 
CCS. 

Between concentric spheres I\ and 5 Dl of BHCS-1, and 
I2 and S 02 of BHCS-2, we reserve overlapping regions, 
which appear shaded in Fig. [TJ These interfaces h and 
S 0i (i = 1 or 2) are not physical boundaries; the bound- 
ary conditions of the fields are not prescribed there. In- 
stead, the value of the field on Ii is calculated from the 
field on BHCS, and the value of the field on S 0i from 
the field on CCS, thus resulting to a smooth potential 
field throughout BHCS and CCS that satisfies the phys- 
ical boundary conditions at the BH boundary and the 
asymptotic region. The significance of the overlap region 
is to decrease the number of iterations to convergence. 
A toy model of this iteration procedure is explained in 
Appendix [DJ 

D. Choices for the Green's function 

In Sec. Ill B[ we discussed a restriction on the choice of 
Green's function due to the KEH iteration using Eq. (fT2")l . 
Any Green's function that meets this restriction can be 
used in the Poisson solver JT2|). The Green's functions 
are expanded in multipoles over the spherical coordinates 
(r, 9, <f>) of each domain. Explicit formulas for the ex- 
pansions are shown in Appendix [Bj In actual numerical 
computations, the summation of multipoles in I is trun- 
cated at a certain finite number L, for which we typically 
choose L ~ 10 in order to resolve the deformation of the 
field $. 

For CCS, we choose the Green's function without 
boundary 

G NB (x,x') = -^—, (16) 
\x — x'\ 

which has the simplest form and picks up the contribu- 
tions from the interfaces I\ and I2. In the volume integral 
of Eq. (fT2"|) over the domain outside of spheres I\ and I2 , 
and inside of S , the function G NB (x,x') is expanded 
in multipoles over the spherical coordinates of CCS. For 



the surface integrals on 1\ and I 2 , G NB (x, x') is expanded 
in multipoles over the spherical coordinates of BHCS-1 
and BHCS-2 respectively. Therefore, the position x cor- 
responding to each grid point of CCS is labeled by the 
spherical coordinates of BHCS, not CCS, in these surface 
integrals. 

For BHCS, when the excision of the computational do- 
main is not used at the black hole, G NB will be chosen. 
However when the computational region is excised in- 
side a sphere Si ± and Si 2 , G NB can not be used for the 
Dirichlct problem. As shown in Appendix [El the I — 
component of V /a G NB (x, x') becomes zero at the inner 
boundary sphere, and hence it can not pick up Dirichlct 
data there during the iteration of Eq. (fT2"j). When the 
black hole boundary condition for a certain field is given 
by Dirichlet data, we choose the Green's function for the 
Dirichlct problem between two concentric spheres G DD 
given in Appendix IB 21 When Neumann data is imposed 
at the black hole boundary, the Green's function without 
boundary G NB may be used. We also coded the Green's 
function between two concentric sphere G ND for which 
the Neumann condition is imposed at the inner bound- 
ary of BHCS Si ± and Si 2 , and the Dirichlet data at the 
outer boundary of BHCS S Ql and S 02 . 



III. METHOD FOR NUMERICAL 
COMPUTINGS 

A. Grid spacing 

Hereafter, coordinate labels (r, 6, 4>) will be used for all 
three spherical coordinate systems, CCS, BHCS-1, and 
BHCS-2 unless otherwise stated. We introduce three 
spheres, S a , Sb, and S c at r = r a , r = H, and r = r c 
respectively, such that r a < r c < r&, in each coordinate 
system. The sphere S a is used as an inner boundary for 
BHCS when the excision boundary is used, which cor- 
responds to S h and S i2 in Fig. [TJ For CCS or BHCS 
without excision, the radius r a = is understood. The 
sphere Sb is the outer boundary of each coordinate sys- 
tem that corresponnds to S a , S 0l , and S 02 in the same 
figure. The sphere S c is located between S a and Sb where 
we change the grid spacing in the radial coordinate. 

The code is constructed to handle non-equidistant grid 
spacing in each coordinate grid. In CCS, the grid starts 
with equidistant spacing from the origin r = to a sphere 
S c , and from there it becomes non-equidistant with an 
ever increasing spacing up to the outer boundary Sb- The 
black hole grids are equidistant from their outer bound- 
aries Sb down to S c , and from that point until the inner 
boundary S a they become non-equidistant with an ever 
decreasing spacing. For the black hole grids, finer grids 
are adequate to have an accurate representation of the 
rapidly varying fields near the hole, while further away, 
where the potentials are changing slowly, larger spacing 
can be used without compromising accuracy. 

We summarize common notations for all three coordi- 
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nate systems as follows; 

r a : Radial coordinate where each grid starts. 

rb : Radial coordinate where each grid finishes. 

r c : Radial coordinate between r a and rj where each 

grid changes from equidistant to non-equidistant 

or vice versa. 

N r : Total number of intervals Ar^ between r a and r^. 
n r : Number of intervals Ar^ between r a and r c . 
n v : Number of overlapping intervals of BHCS to CCS. 
Ng : Total number of intervals A0j for 9 6 [0, 7r]. 

: Total number of intervals A<f>i for <f> G [0, 2ir]. 
d : the separation between centers of BHCS and CCS. 

In particular, we use the following setup for the grid 
spacings. 
Central Grid 



= Aft 
Ar,, 



n r 

fcAr, 



where k > 1. Then we have 



for 1 < i < n r 
for n r < i < N r — 1 

N r -n r + l 



n, 



= Aft 



k-k 



1 



(17) 



Given r a , r^, r c , N r , and n r this is the equation that will 
give us the spacing factor k for the central grid. Note 
that for the central grid, r a = 0. 
Black Hole Grids (I and II) 



Aft 



n 



N r — n r 
Ar t = kAr l+1 

where k < 1. Then we have 



for n r + 1 < i < N r 
for 1 < i < n r 



Aft- 



k - k n - +1 



1 



(18) 



For the angular 8 and <j> spacings for all three coordi- 
nate systems, we usually take equidistant grid spacing. 



B. Finite differencing 

Standard finite difference scheme is applied to eval- 
uate the derivatives of the sources and their numerical 
integrals in Eq. (fT2)) . The derivatives of source terms 
are calculated using fourth order Lagrange formula, and 
the integrals using either trapezoidal rule or fourth order 
Simpson rule in 9 and <j) coordinates and second order 
mid-point rule for the r coodinate. For the surface inte- 
grals at the interfaces S 01 and S 02 in Fig. [TJ the field and 
its derivatives are evaluated from the nearby 64 points of 
CCS as shown in Fig. [5] (a point A on S 0i ), to which the 
fourth order intepolation is applied. 



C. Iteration procedure 

The KEH method at the n-th iteration follows the pro- 
cedure, 




h '4. 



FIG. 2: Schematic figure for computational domain with over- 
lapping grids. 



1) Compute all the source terms in Eq. |T 

2) Call the Poisson Solver (described below) for each 
of the variables and compute their new values 

3) Compare these newly computed values $^™- ) with 
those of previous iteration <£>(™ -1 ). 

• If the difference is less than your accepted error 
=> convergence 

• If not, update $("), according to 



$(") := c $ (n) + (1 - c) 
and go back to step (1). 



(19) 



We conclude convergence of the iteration when the dif- 
ference between two successive iterations becomes small 
as defined by 



2|$(") - $("-!) 



(20) 



where e c = 1CP 8 is taken in typical calculations. The 
iteration usually converges successfully, taking the con- 
vergence factor c to be around 0.5 when a fluid source is 
present. For the binary black hole case, it is even possible 
to achieve convergence with the factor c = 1. 

The Poisson Solver for a potential at the nth iteration 
performs the following sequence of instructions : 

1) Compute the potential and its radial derivative at 
the outer boundary of the black hole grid by inter- 
polating from nearby points of the central grid. 

2) Compute the surface integrals at the outer bound- 
aries of the black hole grids by using the potential 
and its derivative from step (1). 
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3) Compute the surface integrals at the inner bound- 
aries of the black hole grids by using the boundary 
conditions for the potential or its derivative. 

4) Compute the volume integrals inside the black hole 
grids. Add the contributions from steps (2) and 
(3), to obtain the value of the potential inside the 
black hole grids. 

5) Interpolate using points of the black hole grids to 
compute the potential and its derivative on the in- 
terfaces I\ and I2 of the central grid. 

6) Compute the surface integrals on I\ and I2 by using 
results from step (5). 

7) Compute the surface integral at the outer boundary 
of the central grid by using the boundary conditions 
for the potential or its derivative. 

8) Compute the volume integral inside the central 
grid. Add the contributions from steps (6) and (7), 
to obtain the value of the potential inside the cen- 
tral grid. 



IV. CODE TEST 

In this section, we show the results for the convergence 
test of our new code. In the first test we compute the 
Newtonian potential of two spherical masses. Then we 
compute simple models for time symmetric black hole 
data. 

As mentioned in the previous section, the local trunca- 
tion errors of the finite differencing used in our code are 
of order 0(8r 2 ), 0(89 A ), and 0(S(f> 4 ) in each coordinate. 
Also we have a truncation error from multipoles higher 
than L ~ 10. Since the local truncation error at each grid 
point is a linear combination of these, and because of the 
excision used in CCS, we find a non-uniform convergence 
of relative errors as we increase the resolution, as well 
as a non-uniform distribution of the errors in space, as 
shown below. However, overall convergence is faster than 
second order for all cases with a fixed L. 



A. Convergence test for the Newtonian potential 

1. Set up for the test problem 



The Poisson equation (|T0[) with a spherical source of 
the form 



S(r) 



(R 2 - r 2 



if < r < R, 



if r > R, 



(21) 



Type 


Coordinate 


r a 


n 


Tc 


d 


N r 


n r 




No 




L 


SI 


CCS 





100 


3 


— 


80 


40 


— 


20 


80 


10 




BHCS-1 





1.25 





1.5 


30 





6 


10 


40 


5 


S2 


CCS 





100 


3 




160 


80 




40 


160 


10 




BHCS-1 





1.25 





1.5 


60 





12 


20 


80 


5 


S3 


CCS 





100 


3 




320 


160 




80 


320 


10 




BHCS-1 





1.25 





1.5 


120 





24 


40 


160 


5 



TABLE I: Coordinate parameters, and the number of grid 
points for each coordinate system with different resolutions. 
Each resolution is double the one above. The parameters for 
BHCS-2 are identical to those of BHCS-1. L is the highest 
multipole included in the Legendre expansion. 



has the solution 

Rl - 
6 



\R 



3 / r \ 4 1/r 
5 \RJ ~7\R 



8R 3 
105 r 



H0<r<R (22) 
if r > R 



where r is the radial coordinate and R a constant. 

The source (|2"TT) is centered at each BHCS-1 and 
BHCS-2, whose positions in CCS are (r, 9, <fr) = 
(1.5,7r/2,0) and (1.5, tt/2, tt), and the radii of BHCS-1 
(S 0l ) and 2 (S 02 ), extend up to rf, = 1.25. The radii 
of the excised spheres in CCS I\ and I2 are taken as 
r = 1.0. The radial coordinate r of CCS is equidistant 
until r c = 3 and from that point until rf, = 100 is non- 
equidistant, while the BHCSs are equidistant in the ra- 
dial coordinate. The exact potential of two sources is 
a superposition of solutions (f2"2"|) centered at each of the 
two BHCS. For the boundary condition at r = 100 in 
CCS, we put the potential to have its exact value. The 
KEH iteration explained in previous sections is applied 
to calculate the potential $ until convergence is made. 



2. Accuracy of numerical solutions 

In Fig.[3l we show, for two cases, the percentage of the 
relative error between the numerical and exact solution, 



5$ 




[%] := 100 







$0 



rical 



(23) 



In the first case (top panel), the radius of the source R 
is R = 0.5, which is smaller than the boundary radius 
r b = 1.25 of BHCS-1 and BHCS-2. For the second case 
(bottom panel), the source radius R = 1.4 is taken so that 
the sources extend to CCS. The error is plotted along the 
radial coordinate at (0,<fi) — (w/2,0) (see Fig. [1]). Along 
this line (labeled x in the figures) BHCS-1 extends from 
x = 0.25 to 2.75, the source (for R = 1.4) from x = 0.1 
to 2.9, and overlap of CCS and BHCS-1 from x = 0.25 
to 0.5 and from 2.5 to 2.75 in CCS. The errors in the 
interval x £ [0, 4] are shown in the plots. 
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FIG. 3: Percentage of relative error for the source of type 
(|21f) for three different resolutions. Each line from top to 
bottom corresponds to the resolutions SI, S2, and S3 in table 
I. Vertical lines are the location of the boundary of numerical 
domains. Top panel: The source is inside BHCS, R = 0.5 < 
Tb = 1.25. Bottom panel: The source extends outside of the 
BHCS, R=lA>r b = 1.25. 



errors) that we obtain here. 



3. Convergence property of iterations 

In test problems presented in this section, the number 
of iterations required to achieve convergence is about 10- 
20. Generally convergence can be accelerated using (1) 
a good initial guess for starting the iteration, and (2) a 
carefully chosen convergence factor c in Eq. (flT))) . 1 If 
the source is very steep and the convergence factor large 
(near 1.0) the iteration could blow up instead of converge. 
The number of iterations also depends critically on the 
size of the overlap region between the surfaces S 0i and 
Ii (i — 1,2) (see Appendix [D]) • For the results shown in 
Fig. [3J in which the overlap is 20% of the radii of 5*0; , 
the solution is evaluated after 14 iterations with a con- 
vergence factor c = 0.8, starting from a constant (zero) 
potential. For this problem, we can achieve convergence 
with a convergence factor c = 1.0 in 11 iterations with 
the same overlap region. On the other hand, when the 
overlap is set about 2.4% of the radii of S 0i , the iteration 
does not converge with c = 0.8, and it does for c = 0.1 
after 168 iterations. Although the larger overlap region 
is favourable for having the number of iteration smaller, 
the radii of Ii has to be taken large enough to keep the 
number of multipoles used in CCS small. Our choice for 
the radii of S 0i and Ii meets these two requirements. 

We also observed that our method produces the same 
solution irrespective of the value of the convergence fac- 
tor, in the above range 0.1 < c < 1, as long as the itera- 
tion converges. 



B. Convergence test for solutions with excision 
boundaries 



Analytic solutions 



In Fig. [3J the resolution doubles from the top to bot- 
tom (dashed, dotted, and solid) lines in each panel whose 
parameters are shown in Table fl] The errors of the low- 
est resolution, top lines in each panel of Fig. [31 are fairly 
small. For the case with R = 0.5 in the top panel, the 
error drops 1 /4 near the center of the source as we dou- 
ble the resolution, which shows the second order conver- 
gence. For the case of larger source in the bottom panel, 
the volume integration in CCS introduces a truncation 
error that behaves differently from the former case. Re- 
gardless of that the error drops again roughly as 1/4. 

It is remarkable that the potential of such binary 
sources can be accurately computed with the small num- 
ber of multipoles as L = 10 in CCS. In previous work 
[2l| . for binary neutron stars, in which only one domain 
corresponding to CCS was used to compute the field, a 
summation of more than 30 multipoles was required to 
obtain an accuracy of order 0.01% (regarding the relative 



To test our elliptic solver for the case with black hole 
excision boundaries, we consider the following simple so- 
lutions which model one or two black holes. For the met- 
ric 

ds 2 = -a 2 dt 2 + ^fijdx^x 1 , (24) 

where is the flat metric, the hamiltonian constraint 
and the spatial trace of the Einstein equation G a p"f a P — 
0, give 

V 2 yj = and V 2 (cm/0 = . (25) 



When one iterates the fluid variables together with the gravi- 
tational fields, such as a computation for binary neutron star 
equilibrium, the number of iteration may increase as many as a 
few hundreds. 
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These equations have solutions, 



One black hole data without source 



4, 



M 
2~r~ 



and aip = 1 



M 
2r~' 



(26) 



which correspond to the Schwarzschild solution with 
mass M, in isotropic coordinates, V'lr— ><x> — lj an d 

Ot\r — >oo — !• 

We compute the solution (|26[) numerically by imposing 
boundary conditions at the sphere r a — M/2. In order to 
test the code using Dirichlet boundary conditions, we set 
the boundary value at r = r a to the exact value computed 
from (|26p . For testing Neumann boundary condition, we 
take the value of the derivative of ([2"u]) . 

For example, the Neumann condition for aip becomes 



dr 



M 
2^2 



at r = 



M 

Y" 



(27) 



Note that the method of images, [2jJ, is identical to re- 
quiring that ip satisfy the Robin boundary condition 



+ at = 

dr 2r 



at 



M/2 



The boundary condition for the lapse, 



a = 0, 



(28) 



(29) 



yields the solution (|26p . antisymmetric about r = r a . 

One can also construct a two black hole solution that 
satisfies imaging conditions ([28]) and ([29]) at two spheres 
f27l | . Instead of using solutions of this kind, we use the 
Brill-Lindquist solution [28| to Eq. (|25p which is sufficient 
for the purpose of our code test. Writing coordinates of 
BHCS-1 with subscript 1 and BHCS-2 with subscript 2, 
we write a two black hole solution to Eq. ([2"B"]) . 

Mi M 2 Mi M 2 , s 

i> = 1 + + —=- and aip = 1 - - —=-. 30) 
2n 2r 2 ^ 2n 2r 2 V ; 

We set the radii of excision boundaries at 7*1 = M\j2 and 
r 2 = M2/2, and impose either the Dirichlet boundary 
condition, which is given by the values of Eq. ([30]) at the 
boundaries, or the Neumann boundary condition, which 
assigns the values of the derivatives of ip and aip as 



dip 


Mi 


M 2 


<9r 2 


dr\ 


2r\ 


2rl 


9ri 


d(aip) 


Mi 


M 2 


<9r 2 


dr% 


27f 


2rl 


3ri 



at n 



Mi 
~2~ : 



(31) 



at n = (32) 



and 1 «-> 2 for the boundary at r 2 = M 2 /2. The co- 
ordinates ri and r 2 are written in terms of each other 
as 



Type Coordinate 


r a 


Tb 


r c 


d 


N r 


n r 




N e 


iV 


L 


Nl CCS 





100 


3 


— 


80 


40 


— 


20 


80 


10 


BHCS-1 


0.02 


1.25 


0.02 


1.5 


30 





6 


10 


40 


5 


BHCS-2 





1.25 





1.5 


30 





6 


10 


40 


5 


N2 CCS 





100 


3 




160 


80 




40 


160 


10 


BHCS-1 


0.02 


1.25 


0.02 


1.5 


60 





12 


20 


80 


10 


BHCS-2 





1.25 





1.5 


60 





12 


20 


80 


10 


N3 CCS 





100 


3 




320 


160 




80 


320 


10 


BHCS-1 


0.02 


1.25 


0.02 


1.5 


120 





24 


40 


160 


10 


BHCS-2 





1.25 





1.5 


120 





24 


40 


160 


10 



TABLE II: Grid parameters used in one BH test problem. 
The same conventions as in Table [T] are used. 



For each boundary condition, the accuracy of our Pois- 
son solver is examined comparing these solutions to the 
analytic ones. 2 

The Laplace equations ([25]) are solved from the surface 
integrals at the boundaries in our Poisson solver. The 
same equations ([25]) can be rewritten as 



and 



V 2 a 



4, 



.r J <),ri),o. 



(35) 



This form is also used to test the volume integral over 
the source in the Poisson solver. 

In the next two sections Mi/2, M 2 /2 refer to the radii 
of Si 1 and Si 2 of BHCS-1 and BHCS-2 correspondingly. 



2. Convergence test for one black hole solution 

First we treat the problem with no volume sources as 
in equations ([23]) . In Fig. 2] the conformal factor ip is 
plotted along the x-axis. BHCS-1 and 2 are centered on 
the x-axis at x = 1.5 and —1.5 respectively. In BHCS-1, 
we have two surface integrals, one at the inner (excision) 
boundary, , with a radius r a = 0.02 (in conformal 
geometry) and one at the outer boundary, S 0l , at a dis- 
tance rb = 1.25. In the BHCS-2, there is no inner bound- 
ary sphere Si 2 ; we solve for the whole region inside S 02 
without any excised region. In this grid we need only 
to compute one surface integral at S 02 . Finally CCS ex- 
tends to a distance r — 100 and it excludes regions inside 
of two spheres I\ and / 2 which are centered at x = 1.5 
and x = —1.5 respectively and have radius 1.0. The 
overlapping region is one shell centered at x = 1.5 with 
1.0 < r < 1.25 and the corresponding one on the negative 
x-axis. 

In Fig. [5] we show the fractional errors of the conformal 
factor shown in Fig. 3] for three different resolutions in 
table (Til The error of a is of the same order. Since the 
terms calculated using 2nd order finite differencing, such 



- . U2 



2r 2 a sin 6 2 cos + a 2 , (33) 



T2 = \It\ + 2riasin#i cos^i + a?. (34) 



2 In the above two black hole solution, the lapse a takes negative 
value in the neighborhood of each boundary sphere. 
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FIG. 4: Exact and numerical solution for tp on the x-axis 
using the boundary condition ([28} 



„ 10 




FIG. 5: Fractional errors of the conformal factor ip, are plot- 
ted along the x-axis. Lines from top to bottom corresponds 
to the resolutions Nl, N2, and N3 in Table ITTl Vertical lines 
are the boundaries of numerical domains. 



One black hole data with source 


Type Coordinate 


r a 




r c 


d 


N r 


n r 


n v 


N e 




L 


Ml 


CCS 





100 


3 


— 


80 


40 


— 


20 


80 


10 




BHCS-1 


0.02 


1.25 


1.0 


1.5 


120 


112 


8 


10 


40 


5 




BHCS-2 





1.25 





1.5 


30 





6 


10 


40 


5 


M2 


CCS 





100 


3 




160 


80 




40 


160 


10 




BHCS-1 


0.02 


1.25 


1.0 


1.5 


240 


224 


16 


20 


80 


10 




BHCS-2 





1.25 





1.5 


60 





12 


20 


80 


10 


M3 


CCS 





100 


3 




320 


160 




80 


320 


10 




BHCS-1 


0.02 


1.25 


1.0 


1.5 


480 


448 


32 


40 


160 


10 




BHCS-2 





1.25 





1.5 


120 





24 


40 


160 


10 



TABLE III: Grid parameters used in one BH test problem 
with source. The same conventions as in Table U are used. 



CO 




FIG. 6: Fractional errors of the lapse a, are plotted along the 
x-axis. Lines from top to bottom corresponds to the resolu- 
tions Nl, N2, and N3 in Table [Hi] The no-boundary Green's 
function G NB has been used. 



function we get the fractional errors of Fig. [7l With this 
latter choice of the Green's function the error near the 
throat is much smaller. Also we need fewer iterations to 
achieve convergence. 



3. Convergence test for binary black hole solution 



as volume integrals, do not contribute in this solution, 
4th-order convergence can be seen in Fig. [5] as expected. 

Next we solve the same problem but in the form of 
([35)) . Near the inner boundary of BHCS-1, the source for 
the volume integral of the lapse, becomes very steep (ap- 
proximately it goes as r -4 ) and needs more grid points 
than the previous case in order to achive the same order 
of error. For a grid set up as the one shown in Table 
Mil the fractional errors of the lapse are shown in Fig. [5] 
where the no-boundary Green's function G NB has been 
used. The error for the conformal factor is as in Fig. [5] 
since again we don't have any volume sources to inte- 
grate. For the same problem if we use the G DD Green's 



The two black hole solution (|30|) is computed numeri- 
cally solving either set of (|25p or ([3"S")) . In Fig. [51 we plot 
the fractional error, Eq. [53]for the lapse a computed from 
the first set of Eq. ([25)) using the no-boundary Green's 
function G NB and Neumann boundary condition. 

The integral form of the first set ([25)) involves solely the 
surface integrals of Eq. (fT2"]) . Since the surface terms are 
calculated using 4th order finite differences, convergence 
of this order can be seen in Fig. [5] when the grid resolution 
is increased from Tl to T3 in Table HVl Starting from 
a = tp = 1, the solution converges after 15 iterations with 
the convergence factor c = 1.0. Typical CPU time and 
memory to compute the 1/4 of the whole binary black 
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FIG. 7: Same as in Fig. |6] but now G DD is used. 



Two black hole data 


Type Coordinate 


r a 


n 


r c 


d 


N r 


n r 


n v 


N e 




L 


Tl CCS 





100 


2.8 




52 


28 




16 


32 


6 


BHCS-1 


0.1 


1.2 


1.0 


1.4 


32 


30 


2 


16 


32 


10 


T2 CCS 





100 


2.8 




104 


56 




32 


64 


10 


BHCS-1 


0.1 


1.2 


1.0 


1.4 


64 


60 


4 


32 


64 


10 


T3 CCS 





100 


2.8 




208 


112 




64 


128 


10 


BHCS-1 


0.1 


1.2 


1.0 


1.4 


128 


120 


8 


64 


128 


10 



TABLE IV: Grid parameters used in two BH test problem. 
The same conventions as in Table U are used. 



hole region is tabulated in Table Ivl 

In Fig. OH the fractional error (|23|) of the solution to 
Eq. ([55]) is shown for a different choice of the Green's 
functions for BHCS. The numerical integration in radial 
direction, that appears in the volume integrals of (|12|) 
are calculated with 2nd order accurate mid-point rule. 
In all test problems, we found that the largest error ap- 
pears in computing a with Neumann boundary condi- 
tion shown in FigJ5] bottom panel. However, even for 
this case, the error is controlled to give a fractional error 
less than 0.01% everywhere, when the highest resolution 
T3 is used. We have tested different combinations of the 
Green's functions with boundary conditions, and found 



Type CPU time/iteration [s] 


Memory [MB] 


Tl 0.14 


25 


T2 1.0 


73 


T3 8.6 


236 



TABLE V: Typical CPU time and memory used for the BH 
calculations. The use of equatorial and ir rotation symmetries 
reduce the number of grid points indicated in Table HVl by a 
factor of 4. Note that the computational costs approximately 
scale linearly with respect to the total number of grid points 
which is 2 3 = 8 times at each level T1-T3. Opteron 2 GHz 
with Portland fortran compiler is used. 
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FIG. 8: Fractional error of the lapse is plotted along the x- 
axis. Lines from top to bottom correspond to the resolutions 
Tl, T2, and T3 in Table IrVl A set of Eq. fl25j is solved using 
the no-boundary Green's function and by imposing Neumann 
boundary conditions. Vertical lines correspond to the bound- 
aries of the numerical domains. 
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9: Same as Fig. [S] but for Eq. (gSJ. For the BHCS, 
D is used in the top panel and G ND in the bottom panel. 
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x 

FIG. 10: Fractional errors of the lapse a are plotted along 
the x-axis for different values of the highest multipole in the 
Green's function. A set of Eq. (|35fl is solved using the Green's 
function G ND with resolution T3 in Table |TV] 

similar or better convergence results. (The fractional er- 
ror of the conformal factor ip, on the other hand, scales 
in 4th order since the equation for ip does not involve the 
volume integrals. ) 

Finally, the convergence test for the number of mul- 
tipoles summed in the Green's function, L, with a fixed 
resolution (T3 in Table llV|) . is shown in FigdOl Changing 
L from 4 to 12, convergence is achieved around L = 10. 
This is a dramatically small number compared to the 
binary neutron star case [2l[ , for which L = 30 ~ 40 
multipoles are required because only one domain (cor- 
responding to CCS here) was used for computing the 
gravitational fields. 



V. AN EXAMPLE FOR BINARY BLACK HOLE 
INITIAL DATA 

To conclude the test for our new numerical code, we 
show an example of binary black hole initial data; a bi- 
nary black hole solution of IWM formulation. This is to 
demonstrate that our new code can produce the binary 
black hole initial data with nonzero angular momentum, 
and locate the apparent horizon using the method de- 
scribed in Appendix [X] at the same time. Writing the 
spatially conformally flat metric 

ds 2 = -a 2 dt 2 + tp i f lJ {dx l + f3 l dt)(dx3 + ftdt), (36) 

in a chart {t,x 1 }, the 5 metric potentials, the conformal 
factor ip, the shift (3 a and the lapse a are solved from the 
Hamiltonian constraint, momentum constraint and the 
spatial trace of the Einstein's equation, respectively. As 
shown in Appendix [Cl all these eqations are written in 
elliptic form. 

At the black hole excision boundary, certain Dirich- 
let data is imposed to ensure that the apparent horizon 



Type r a i/jb «b fi Ob 
Bl 0.1 3.0 1.0 0.3 0.0 

TABLE VI: Parameters for the boudnary conditions. Except 
for the value of r a , the parameter set of T3 in Tabic HVl is used 
for the computation. 

appears outside of the boundary sphere. For simplic- 
ity, we do not intend to impose certain physically moti- 
vated boundary conditions (see below) . Dirichlct bound- 
ary conditions are given to all variables {ip, a, f3 a } as 

ip = ips = constant, (37) 
a = as = constant, (38) 

/3 a = -n<p a c -n B <p%, (39) 

at the excision sphere r = r a of BHCS. For the boundary 
value of the conformal factor ife, a constant is chosen 
large enough to form apparent horizons near the excision 
spheres. For the lapse qb, we also assign a constant value. 
The boundary condition for the shift vector assigns a 
momentum and a spin to each hole. Here, 4>q and <f>^ are 
the basis of <f> coordinate of CCS and BHCS, respectively. 

In Figs. [Til and [T2l we show a solution with the bound- 
ary parameters shown in Table IVIl All potentials are 
smoothly joined across the overlap of CCS and BHCS. 
In Fig[T21 thick dotted circles right outside of the excised 
sphere (thin black circles) are the apparent horizons lo- 
cated using the method in the Appendix [XJ 

VI. DISCUSSION 

The KEH iteration using the Green's formula Eq. (JX3J) 
is applicable to solve various types of partial differen- 
tial equations with nonlinear sources. For example, this 
method has been applied to solve the helically symmetric 
scalar field and binary neutron stars 29]. In this work, 
the equation for the scalar field is written in the form 
of Helmholtz equation, and the half-advanced + half- 
retarded Green's function is used in Eq. (fT2")l to compute 
a standing wave solution iteratively. We plan to com- 
pute the helically symmetric binary black hole/neutron 
star solution using the coordinate systems and iteration 
scheme presented in this paper. 

Black hole singularities on an initial hypersurface are 
avoided either by excising a numerical domain in the 
neighbourhood of the singularities, or by using punctures. 
In the previous Sec.fVl we applied rather crude boundary 
conditions at the black hole excision sphere. For physi- 
cally motivated excision conditions, one imposes appar- 
ent horizon or isolated horizon boundary condition to the 
conformal factor at the excised sphere so as these spheres 
become automatically horizons [3(J [M[. Alternatively 
one can use the puncture method to produce accurate 
initial data as the ones used in binary black holes evolu- 
tions 0, [f| [EJ3, HI, and black hole/neutron star binary 
simulations [321 ]. 
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FIG. 11: Plots for a (top) and [3 V (bottom) along the x-axis 
which pass through the excised region of the two holes. Ver- 
tical lines are the boundaries of each computational domain. 
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APPENDIX A: APPARENT HORIZON SOLVER 

A method for locating the apparent horizon is de- 
scribed in this Appendix. Our method is a modification 



FIG. 12: Contours of the lapse a (top) and the vector field 
of the shift f3 a (bottom) for the same model as Fig llll Thick 
dotted circle in each figure is the apparent horizon. 



of [331 ] and it will be one of the simplest apparent hori- 
zon finders without any symmetry restriction. For other 
works on apparent horizon finders, see e.g. [34], |35| and 
references therein. 

An apparent horizon A is defined as the boundary of 
all trapped regions on a spacelike hypersurface E, where 
the expansion ■& of the outgoing null congruence £ a or- 
thogonal to A vanishes. Introducing a null foliation 7i u 
whose normal is £ a , where u labels a family of null hy- 
persurfaces, the apparent horizon A is the intersection 
A = Ho HE. To each two dimensional surface TL U (~l E, 
the spatial unit normal vector s a is associated, and the 
outgoing null vector field £ a can be written with a func- 
tion /, l a = f(n a + s a ), where n a is a timelike normal 
to E. Then, writing a projection tensor onto a surface 
orthogonal to s° as e a b = 7 a h — s a Sb, the expansion 1? is 
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written 

= e^VJp = f(D a s a - e a(i K afi ) (Al) 

Hence, the equation 

D a s a - e ab K ab = 0, (A2) 

is satisfied at the apparent horizon A. 

We introduce a family of level surfaces parametrized 
by a spherical coordinate surrounding a black hole in the 
form , 



the equation for the apparent horizon (|A2|) is rewritten 



F:=r-R h {6,cj>), 



(A3) 



where F — coincides with the apparent horizon. 3 The 
spatial normal s a is proportional to the gradient of F, 



_ DgF 
S " \\DF\\' 

where the norm \ \DF || is defined by 

||£>F|| = Vl ab D a FD b F 



(A4) 



(A5) 



We rewrite Eq. (|A2|) as an elliptic equation for the level 
surfaces in the conformally related geometry. Intro- 
ducing quantities weighted with the conformal factor, 
e a b ■= il> A e ab , K ab := ip^Kab, s a := ip 2 s a , the norm (|A3)) 
is transformed 



and 



DF\ \ = v/r 2 Jj ab D a FD b F = i{j- 2 \\DF\ 



D a F 



\DF\ 



(A6) 



(A7) 



Multiplying Eq. (|A"2)) by the factor tf; 2 \\DF\\, we obtain 
an elliptic equation for the apparent horizon, namely 



AF + S = 0, 
S :=D a In- 



(A8) 



V> 4 



-D a F-^\\DF\\K ab e ab , (A9) 



\DF\ 



which is satisfied on the surface r — Rh(6, 4>)- Separating 

o ~ 

the Laplacian associated with the flat metric, A, from A 
associated with the conformal metric, ^ ab , we have 

AF = AF + h ab D a D b F - j ab C r ab D c F, (A10) 

where the second and third terms have been moved to 
the right hand side. Defining 



AF 

O 

A H 



- -\k H R h = --^ ( A H - 2)R h ,{Ml) 
r r l Rr 



sm9d9 { 86 



1 d 2 
sin 2 6 86 



(A12) 



3 The level surface F may not coincide with the intersection WuHS, 
except at A. 



(A H -2)R h = S, 
S := Rl(S + S), 



(A13) 
(A14) 



S := -h ab D a D b F + r b C c ab D c F (A15) 

Terms in S vanish for spatially conformal flat geometry. 
We find a solution to Eq. (|A13j) 



Rh(x) 



4tt 



d 2 xG(x,x')S(x'), (A16) 



H 



where the coordinates x here represents (9, 0), and the 
function G(x, x') is given in terms of Legendre expansion, 

G(xx >) V 2£+1 f E (Z-rn)l 

x P; n (cos 6)Pl n (cos 6') cos m(<f> ~ <j>'). (A17) 

The same discretization as in BHCS, and an iteration 
similar to KEH method are applied to Eq. (|A16[) (sec 
Scc lIIII) . The 4th order Lagrange formula is used for fi- 
nite differencing the source and for numerical integration. 
The iteration converges typically in 30 iterations, whose 
CPU time is negligible in the computation of initial data. 



APPENDIX B: MULTIPOLE EXPANSION OF 
THE GREEN'S FUNCTIONS OF THE 
LAPLACIAN. 



In our Poisson solver (|12[) one must choose appropriate 
Green's functions to meet the boundary conditions im- 
posed on each of the field variables. In this Appendix, we 
present explicit forms of those used in the preceding sec- 
tions. They are the Green's function without boundaries, 
G NB (x, x'), and two Green's functions with boundaries 
on two concentric spheres S a and S b at radius r = r a 
and r = r bl where r a < r b ; one of them imposes Dirich- 
let conditions on both S a and S b , G BB (x,x'), and the 
other imposes Neumann condition on S a and Dirichlct 
condition on S b , G NB (x,x'). All Green's functions, rep- 
resenting G(x,x'), are expanded in multipolcs, in terms 
of the associated Legendre functions P™ (cos 6) in spher- 
ical coordinates (r, 6, 4>), 



G{x,x') = ^2g e {r,r')^2 e, 



(l-m)\ 



(£ + m)\ 

x pp (cos 6 ) P} n (cos 6' ) cos m ( - <j>' ) , (B 1 ) 
where the coefficient e m is defined by 



1, for m = 0, 

2, for m > 1, 



(B2) 



and hence the difference appears in the radial part of the 
Green's function gi(r,r'), [36| . 
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1. Green's function without boundary G NB (x, x') 

For the Green's function of the Laplacian without 
boundary G NB (x,x'), the radial part gf B (r, r'), is de- 
fined by 



,NB 



(B3) 



where r> :— sup{r, r'}, and r < :— inf{r, r'}. 

When the Green's function G NB is applied to BHCS 
with two concentric boundary spheres S a and Sb located 
at r — r a and rb, respectively, the following values are 
used to evaluate surface integrals of Eq. (fT2]) , 



9 f B (r,r a ) 



9? B (r,r b ) = 



^+1 ' 



e+i ■ 



and, since V' a G NB (x, x')dS' a = +<9 r 'G r dr'dfl' 



dr'9f B (r,n) = -(*+!)- 



,1+2 ' 



(B4) 
(B5) 

(B6) 
(B7) 



Note that the form of (|B6|) indicates that the Green's 
function G NB does not pick up £ = mode of the Dirich- 
let data at the sphere S a (r = r a ). 



2. Green's function for the region between two 
concentric spheres with Dirichlet conditions, 

G DD {x,x') 

The Green's function G UD (x,x') is a solution of Eq. 
(fl~Tj) in a region between two concentric spheres S a and 
Sb with radius r = r a and r = Tb [r a < rb) where Dirich- 
let conditions are imposed. Its radial part gf u (r, r') as- 
sociated with the £th mode is written 



n 
t+i 



21+1 



-i -1 



^+1 



£+1 



(B8) 



By construction <7 BD (r, r') vanishes on the two spheres 
S a and S b , 



0. 



5 BD (r,r b ) = 0. 



(B9) 
(BIO) 



The derivatives that are used to compute the surface in- 
tegral in Eq. (fl"2j) . are 



dr>9? u (r,r a ) 



21+1 



x(2£+l) r -j^ 



(?) 



d r >g? D (r,n) 



l - 



2i+l 



,(B11) 



x(2£+l)^ 



e+2 



(t) 



t+i 



(B12) 



at S a and Sb, respectively. 



3. Green's function for the region between two 
concentric spheres with Neumann and Dirichlet 
conditions, G ND (a;,a;') 

Similarly, G ND (x,x'), is the Green's function between 
S a and Sb, where Neumann data are imposed on S a 
and Dirichlet data are imposed on Sb- Its radial part 
gf D (r,r r ) associated with the £th mode can be written 



i 



i 



£+l\r b 



21+1 



^+1 



£+l\r < 
r> 

n 



t+1 



(B13) 



The values of g^ at the surfaces S a and Sb become 

/ .. \ - 1 i 1 

JMD/ 



9 r(r,r a ) 



1 + 



2£+l r l a 
X £ + 1 r e b +1 

g? D (r,r b ) = 0, 

and its radial derivatives, 

d r >gf D (r,r a ) = 

dr'9? D (r,r b ) = - 



+ 1 \n 



2£+l' 



x(2£+l)^ 



£+2 



£ /r..\ 



t(t) 



(B14) 
(B15) 

(B16) 
(B17) 



The values of d r >g^ (r,r a ) and gf D (r, rb) vanish by 
construction. 



APPENDIX C: 3+1 DECOMPOSITION FOR THE 
EINSTEIN EQUATIONS 

In the usual (3+1) decomposition of the Einstein equa- 
tions the spacetime metric is written as 



ds 2 = g a pdx a dx^ 



15 



a dt + 7i,(eb I + f3 l dt)(dx 3 + (3 3 dt) (CI) the Ricci tensor becomes 



where g a p, 7y, are the 4D and 3D metrics, while a and 
/3\ are the lapse scalar and the shift vector respectively. 
The Riemannian 3-metric 7^ on a hypersurface E is iden- 
tified by the 4-tensor 



lap = g a p + n a np 



(C2) 



where n a — — a\7 a t is the unit future pointing normal to 
the hypersurface S. The indices of j a p can be raised ei- 
ther by 7 Q/3 or by the full metric g a @ and that 7"^ projects 
vectors onto the subspace orthogonal to n a . Note that 
g a p and "/ a p differ only on the time-time component while 
g a P and 7°^ have identical only the space-space compo- 
nents ( ff Q % 7 = 5% 7 « 7j . fc = 6^, 7 Q,3 7 / 3 7 ¥> The 
covariant components of the shift are (3j — "fij(3 l and the 
components on the normal vector are 

n a = (a, 0, 0, 0) and n a = (-,—) (C3) 

a a 

The extrinsic curvature is 

K a p = -D a np = -7 Q Q 7/ f Va'tip' = - ^£ n J a p (C4) 

where V, D are the covariant derivatives associated with 
g a p and i a p respectively. The Einstein equations can 
now be split into the constraint equations 



Ru = A(A t 2) DitpDjip - A( ^ . 2 2) jijD m tpD m ^ 



4*0 2 
A 



4V> 2 

(DiDjifi + ^ijD m b m ^) + ilij (Cll) 



2tp 

and the Ricci scalar 
^^^1^^^"^^' (C12) 



Now the Hamiltonian (|C5[) and the momentum constraint 
(IC6I) can be written as 



a , i> ^ A — 4 ~ 
Atp- ^-Tl 1 



A+l 



(C13) 



2A-< 4^ AW+^r-^-^ 

2A~^ = Y~ 

+ — 1\ '■' D , <. • - ^KfW^ 

-^DjK = 8nf . (C14) 



-Dfc is the covariant derivative with respect to 7^ and 
A 



A = D l Di. Note also that Dnjj = DiLOj — C™jU> m and 



R /\,,/\'' • l\' = 167T/9 

Dj{K ij -f j K) = 8nf 



and the evolution equations 



d 7y 
dt 



= -2aK i:j + 2Di(3j + 20,(5, 



(C5) 
(C6) 



(C7) 



= allij - DiDja + a{KK l3 - 2K im K j 

+K mi D 3 f3 m + K m jDi(3 m + rD m K l3 (C8) 
1 



hra[T ij + - r f ij (p-T™) 



where p = T a pn a n^ and j a = —Tp^n^j" 13 are the en- 
ergy and momentum density respectively as seen by an 
observer with four velocity n a while IZtj , 1Z, are the three 
dimensional Ricci tensor and Ricci scalar on the hyper- 
surface S. From the two evolution equations we can find 
the time derivative of the trace of the extrinsic curvature 

d t K = aK-Aa + aK 2 + f3 l D i K - 8^7^ (C9) 

where Pij = Tij + ^"fijip — T^) are the source terms and 
A = D i D i . 

With a conformal transformation of the form 



7u = ^ A lij 



(CIO) 



1. Field equations for the initial data 

Since we are searching for quasi-equilibrium states we 
assume the existence of a Killing vector 

e = e + c=(0 + n{^y = {lxi) ( oi5) 

where = £l(d / 'dip) 1 , and is a constant representing 
the orbital angular velocity. In the presence of £ Q the 
spatial metric 7^ and the extrinsic curvature Kij satisfy 

£au = and £ $ K v = • ( C16 ) 

From the first equation of (|C16|) we get 



dt-Jij + ACj + DjQ = 



(C17) 



thus the time derivative of the spatial metric is associated 
with the spatial derivative of the rotational vector £*. 
With the help of the evolution equation of 7 y , (|C7|) we 
can obtain an expression for the extrinsic curvature by 
eliminating the time derivative of the three metric and 
therefore cast the initial value equations in a form that 
has no time derivatives. We find 

Ka = ^-(D.uji + DiUi) and K = -fl.w' (C18) 
2a a 
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where u>i = [3i + Q is the comoving shift. Since 



1 



KijKV = — (Lw)^ (La;)« + -K 2 (C19) 



where 



(Lw)ti = Awj + Djuji - - ltl D m Lo m (C20) 



In the momentum constraint, the last two expressions, 
equations (|C28p and (|C29[) , come from the fact that D^K 
involves the second derivative of the comoving shift uoi as 
follows 

DiDjU) j = aDiK + KDia - A f -^oV A/^J (C30) 



and (Lw)y(La;) y = (IiD) tf ( La))^ the hamiltonian ([C13|) 
and momentum constraints (|C14|) are written 



ib - A - 4 ~ - ■ ?/> A+:L - ~ • ■ 



2A Atb 
?A A+1 ^ 2 87r j oV A+1 



8Aa 2 



A 



(C21) 



4a -, 



o 

Also by using the fact that 



A - 



Aa = ib~* ( Aa + —D i ipD' l a 



(C22) 



(C23) 



we can rewrite equation (|C9[) in the conformal geometry 
as 

Aa = —^D i ipD i a + ip x ( ^(La%(La)) ij 
2^ \4a 



c^A-fr + 47ra(p + T) (C24) 



aK 2 



For the binary black hole case, the sources p, Tij and j z 
vanish, and take A = 4. Under these assumptions, our 
system of equations is 



Aa = ^-(Ld»)y(LJ)) y + + a) l A# 



-Diib&a 

ib 



(C26) 



Aw, = —DiDjUp - 7t y AV + A In \JjL ) (La; ), , 



4a - 

DiK 

3 

= -Kijtf +D< In ( ^ ) U^);, 



(C27) 



2 



if ^ 



+A ( ^ D ^) - yAa + «A# (C28) 

-7?y^'+j;'lii(^j (La)),-, 
4if - 

-A ( ^ D i^) ~ — D * a + AAw 3 (C29) 



APPENDIX D: TOY MODEL FOR 
IMPROVEMENT OF CONVERGENCE BY 
OVERLAP REGION 

To analyze the improvement of the rate of convergence 
achieved by the overlap region, we consider a simple 
model to calculate the potential of a point mass M us- 
ing two overlapping concentric spherical grids as shown 
in Fig. Q21 bottom panel. The first grid extends from 
the coordinate center (r = 0) to the surface Si, and the 
second grid from S2 to infinity (or practically to a large 
distance). In Fig. 1131 top panel, there is no overlapping 
region, Si = S2 =: S. The potential in region II <I>n, can 
be calculated from the surface integral of the interface 
5*2, as 



M + m 



(Dl) 



where m is the error inherited from the initial guess for 
<I> given at the boundary S 2 ■ The potential of the region 
I $1, is a sum of a volume and a surface integral at Si, 



M 

$1 = he. 

r 



(D2) 



Again e is the error inherited from the initial guess for 
$ at Si. From continuity of potentials at Si, $i(i?i) = 
$ii(i?i), therefore 



$1 = 



M m 



(D3) 
(D4) 



The next step of iteration is to use this value to fix $n- 
In region II, the potential at R 2 must satisfy <l>ii(i?2) = 
<I>i(i?2), whence 



M - 



Iterating this procedure n times, we will have 



M {n) = M + : 



Ri 



(D5) 



(D6) 



namely, the error after the n-th iteration is m{R2/R\) n - 
Analogously for the £-th multipole component, writing 
the actual value as A, and the error in region II as a, 



A 



r e+i 



(D7) 
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The solution to region I will be 

A 



Br" 



(D8) 



where the second term is again an error. The boundary 
condition at R\, $i(i?i) = $n(i?i) yields 



B = 



a 



Ri 



2£+l ' 



(D9) 



and from continuity at i? 2 , $11(^2) = $1(^2), 

, 2£+l 

$„ = ^ . (D10) 



Hence the higher multipole converges faster after the n-th 
iteration, 



aw =a + , 



R2 



i(2£+l) 



(Dll) 



The boundaries Si and S2 are used to communicate 
the information of the physical boundary conditions im- 
posed at the asymptotic region, and the inner excision 
boundary, as well as the source from each region to the 
other. At the n-th step of iteration, the values of <&(x') 
and d$(x')/dr' in the surface integrals in Eq. (jXSJ) are 
calculated from the potential of (n — 1) step of the iter- 
ation. 

It is possible to achieve communication between the 
two regions without the overlap region, by mixing the 
values of the potentials <&i and $11 of the two regions. 
This is done when we calculate the value of d$>(x')/dr' 
at the boundary as shown in the top panel of Fig. fT3| 
choosing say the values at grids £3, xi, y±, 1/2, TJz- In this 
way convergence to a correct solution is again obtained, 
but the number of iterations increases approximately ten 
times, even for simple toy problems presented in Sec lIVI 

In the actual binary calculation shown in Fig. [21 when 
the values of <& and d<&{x')/dr' on S 0l , are intepolated 
from CCS, some of these points do not belong to the 
computational domain of CCS. The smaller the overlap- 
ping region the more of these points exist. (For example 
at the point A in Fig. [21 when we interpolate the nearby 
CCS points we find that the point in the lower left corner 
does not belong to CCS.) In such case we interpolate the 
nearby BHCS points to find the value of the potential. 



APPENDIX E: CONVERGENCE OF THE 
ITERATION 



In section III Bl we have said that a requisit of the 
Green's function for the KEH iteration to achieve con- 
vergence is that all multipole components of V a G(x,x') 
should not vanish at the boundary for solving Dirichlet 
problem, and similary G{x,x') should not vanish for Neu- 
mann boundary conditions. In this Appendix we use the 
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FIG. 13: Top, solving the poisson equation on two non over- 
lapping grids. Bottom, solving the poisson equation by using 
two overlapping grids. 



fixed point theorem to illustrate the convergence (or not) 
of the KEH method for simple spherically symmetric case 
with the flat Green's function. 

Let's start with the following problem: 

V*ip = , r > r a with ^ + = atr = r a , (El) 

or 2r 

and linir— too <f(r) = 1, whose solution is (p = 1 + r a /r. If 
we consider the map defined by 



G{x,x)—- V {x) — 



dS' (E2) 



on a suitable Sobolev space, where S a is the surface r 
r a , we can show the following: 
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claim: If G(x, x') — \ x l x ,\ and the functions ip(r) satisfy 
— - on 5„ then $ has a fixed point. 



dip | 
dr ~ 2r 



Indeed if ip\ and (f2 satisfy the above boundary condition 
so is ifi — if2 therefore 



d($Oi),i>(vj 2 )) = sup{|$((pi) - &((p 2 )\} 



d(ipi — 1P2) 



(E3) 



where for the term with g~> we used the boundary 
condition and for the term with <pi — <p 2 we used the 
fact that it depends only on r thus can be pulled out 
of the integral which when calculated gives zero. There- 
fore <I> has a fixed point which can be found if we take 
an initial value ipa and then compute the sequence {</? n } 
with (fin+i — ^{fn) (KEH method). By doing so and 



adding the contibution from infinity which is 1.0 we get 



Vn+i{r) = 1 + 



r a tp n (r a ) 
2r 



which tends to 1 



Now if we change the boundary condition to (p = at 
r — r a the solution turns out to be tp — 1 — r a /r. In this 
case the above argument breaks down and as we have 
seen the KEH iteration fails. The above argument gives 
us 



d($(<pi),$(<p 2 )) < r a 



d(ipi - ip 2 ) 



dr 



(E4) 



and nothing guaranties that $ will have a fixed point 
any more. Actually tp n +i(r) = 1 - r -f \^§f) thus 

V / r—r a 

starting from any constant value, the sequence is stuck 
at 1.0 and this explains why our code gives everywhere 
the value 1.0 with the Dirichlet boundary condition. 
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